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Abstract 

Interface conditions for coupling the domains in a physically motivated domain 
decomposition method are discussed. The domain decomposition is based on an 
asymptotic-induced method for the numerical solution of hyperbolic conservation laws 
with small viscosity. The method consists of multiple stages. The first stage is to 
obtain a first approximation using a first-order method, such as the Godunov scheme. 
Subsequent stages of the method involve solving internal-layer problems via a domain 
decomposition. The method is derived and justified via singular perturbation tech- 
niques. 


1 Introduction 

This is a report on a preliminary investigation of conditions for the interfaces between sub- 
domains when solving partial differential equations. The analysis for the method is a combi- 
nation of asymptotics and numerical analysis. The result is a physically motivated domain 
decomposition method where different partial differential equations may be solved in different 
domains. Since different modeling equations are in different subdomains for the same prob- 
lem, we call this heterogeneous domain decomposition. The numerical treatment of interface 
conditions between the subdomains must be addressed. The approach here is to examine 
the physics reflected in the numerical method used within the subdomains and guarantee 
that this same physics is reflected in the interface treatment. 

The method is best suited to partial differential equations that contain regions of singular 
behavior. A typical situation is when there are narrow regions where the variation in the 
solution is large. Such regions are called boundary layers or transition layers depending on 
whether they are near a boundary or inside the interior of the domain. Examples of such 
situations are laminar flow of a slightly viscous fluid or combustion with high activation 
energy. Classical schemes applied to these types of situations generally fail to correctly 
describe the behavior inside the layers. This difficulty is overcome by utilizing asymptotic 
analysis that reflects the physics of the problem. Here we present and motivate the domain 
decomposition method, but the details of the analysis are presented elsewhere [7]. 

There have been some intersting results regarding interface conditions for heterogeneous 
domain decomposition where Euler equations are coupled with Navier-Stokes equations [9], 
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and where viscous and inviscid equations where coupled [2, 4]. Many of the basic ideas 
relating to asymptotic analysis and numerical methods that utilize domain decomposition are 
found in [10]. These ideas were incorporated into a parallel numerical method in [5]. Specific 
application to conservation laws have been developed in [1]. There are other important" works 
in these areas-these references are only a small sample of the literature. 

The coupling of the problems in the subdomains is based on a balance of the flux across 
the interface. Each subdomain is treated as a control volume, and the flux into and out-of 
the control volume is balanced. This is similar to the flux-differencing methods used within 
the subdomains. The result is a numerical method with no visual artifacts. This numerical 
treatment of the interface is an extension (to heterogeneous domain decomposition) of the 
work by Osher and Saunders [11]. We expect extension of this method for the interfaces to 
work for two dimensional heterogeneous domain decomposition, since it was used for a two- 
dimensional homogeneous domain decomposition method that utilizes adaptive refinement 
[ 8 ]. 

2 Problem Setting and Domain Decomposition Mo- 
tivation 


Consider the Cauchy problem 


f + inu) 


*7(2,0) = V{x ) 


' dx 


{ p ( u )-i£) for 

for x G IR. 


( 2 . 1 ) 


Here the solution U G fll n is a vector- valued function with n components, the domain is 
0 = IRx]0,T[ and e << 1 is a small parameter. 

We assume that V is piecewise smooth. We also assume F and P are regular functions 
of U. We suppose that P is a suitable viscosity matrix [3] for the shocks of the associated 
inviscid problem 

s sr + &F(V°) = o for (M)efl 

( 2 . 2 ) 

k *7°(a:,0) = V(a:) for x £ 1R. 

Namely, a shock-wave solution to (2.2) can be obtained as a limit of progressive wave solutions 
of (2.1). Problem (2.1) is a parabolic-hyperbolic singular perturbation problem driven by 
( 2 . 2 ). 

The regions where the solutions to the associated inviscid problem fail to be good ap- 
proximations to the solution of the full problem are the regions where we use a subdomain 
to localize the behavior of the solution. Thus, we have two types of domains. The first type 
of domain is located where the regular expansion v 


pouter = ^0 + cU l + C 2 V 2 + ^ . 


(2.3) 
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for U is valid and the solution is smooth. The second type of domain is where the solution 
exhibits singular behavior and the regular expansion for U is no longer valid. 

We substitute {7° s uter in the differential equation of (2.1) and use identification in c to 
obtain that U° must be a solution of (2.2). The inviscid problem (2.2) has many weak 
solutions; it is possible to uniquely define U° by considering the problem that governs U 1 

[?]• 

The failure of the regular expansion is reflected by some of the terms in the PDE governing 
U° being significantly larger than other terms. Typically, the term RHS(U°) will become 
unbounded as the small parameter e tends to zero. For finite e, a large RHS(U°) would 
indicate that the region should be covered by a subdomain in which we apply techniques 
designed to capture the singular behavior of the solution. We describe how to use a measure 
of the numerical approximation of RH S(U°) to place the subdomain boundaries in a later 
section of this manuscript. 

; # 

2.1 Problem in the Singular Region 

So that we can handle the regions where solutions to problem (2.1) contain shocks that 
interact with other singularities we use a brute force approach that will capture all possible 
behavior of the solution. The approach is to use the coordinate system 



in the regions with shocks. We will present and motivate the domain decomposition method, 
but the details of the analysis are presented elsewhere [7, 6]. Under this transformation the 
PDE that governs the solution becomes 

(2 - 4) 

where f/(£,r) = U(x,t). This is the equation that is solved in the singular region. 

This scaling is most appropriate for regions where shock-layers are interacting with other 
non-smooth physical phenomena. Because the transformation a priori resolves all of the 
physics. This is reflected by all of the terms in (2.4) having magnitude of order unity or 
smaller. In general, this method is overkill, similar to using a shotgun to dispatch a housefly. 
We choose to study only this brute-force approach so that we concentrate on one type of 
interface. Other treatments that include more of the physics are possible [7]. They can result 
in more efficient numerical methods than the one discussed here. 

The boundary condition at the interface is to impose that the viscous equation from 
problem (2.1) be the model at the interface between the subdomains. The computational 
implications of this condition is discussed in §4. 
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3 Conservative Discretizations 


It is important for the discretization techniques to satisfy a discrete conservation relation. 
One can verify that if the discretizations can be written in the form 

z i +l = z i “ M^t+i/2 “ 1 / 2 )? 

then the method satisfies the appropriate conservation relations. Here we use flux differencing 
methods based on a finite-volume formulation of the problem. 

We will discuss the differencing method for the outer region subdomain where the solution 
is smooth first. Let Wq be the discrete numerical approximation to U°. We use a first-order 
finite-volume method. This method assumes that the value W q t * is an approximation to the 
average of the desired function U° over the spatial interval ]zi_i/ 2 , ®;+i/ 2 ] a t time t = &At. 
The method can also be categorized as a flux differencing technique since the general form 
of the discrete analogue to the original PDE can be written 

Wtf/ 1 = Ki - W+1/, - lt,„) ■ (3.5) 

where 

•*?«/, *F( Kn/i)- 

Here the fluxes are based on the first-order Godunov scheme; thus, the flux fj for com- 
ponent wj of Wq is approximated as 


fU 1/2 = \ [fi(K.) + - a.W, + . - WW] 


(3.6) 


where a* is an approximation of the upper bound on the local speed of sound. 

The discretization that is used for the numerical method in the shock-layer region is 
a modification of the treatment used for the outer region. We have used a coordinate 
transformation that creates a smooth problem for this subdomain. Let Wq be the first 
order numerical approximation to U. Let W^ be an approximation to the the average of 

the desired function U over the spatial interval ]£;_ i / 2 i £ i + i / 2 ] at ^ me r = fcAr. The flux 
differencing technique is 

(3.7) 


Wf" = Wt- A (f? +1/! - Ft,,,) 


frit 
+ 1/2 


where 


frj< 

1 1+1/2 


nHuJ 


- 1/2 ) 


l/2> Qf 

The particular discrete form for each component of the flux is obtained using a formula 
similar to that of Equation (3.6). 

We are not restricted to this particular numerical discretization; however, the numerical 
treatment of the interface will possibly need to be modified for different numerical treatments 
of the problems within the subdomains. 

One can verify that the flux differencing methods given above satisfy the discrete con- 
servation relation. What remains is to formulate the conditions at the interface so that the 
relation will be satisfied globally. 
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4 Treatment of the Interface 


Using the shock-layer coordinates with A£ = C Ax will result in G/e points in the shock- 
layer for each point in the outer region. Here, a typical value for t is .01; hence, this results 
in a radical grid refinement for the shock-layer. For the numerical method, since there will 
be many grid points in the shock-layer for each point in the outer-region, we will refer to 
the shock-layer grid as the refined grid, and the outer region grid will be called the coarse 
grid. The temporal coordinate will also be stretched, resulting in the situation outlined in 
Fig. 4.1. 


x i- 1/2 


x i+l/2 


time 



6/2 


space 

Figure 4.1: Interface at the left boundary 


4.1 Flux Treatment of Interface 

As in [11], we view the interface treatment as a predictor-corrector method on the coarse 
mesh. We start at time t = t k . The coarse-grid values are defined everywhere, and are the 
average of the corresponding fine-grid values when the coarse-grid volume element is within 
the fine-grid region. 

The steps for the first order method are outlined in Algorithm 1 below. At time step 
k, the shock-layer has N(k ) points in the interior of the region and a ghost point on each 
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For k = 1, 

I. March Wo from 4-i to 4 based on scheme (3.5). 

II. Detection. 

A. Compute the residual on the coarse mesh. 

B. Mark regions that should be refined. (Let this be the region between x% L -i / 2 
and x iR+ 1 / 2 . 

A. Modify shape of refined region. 

III. March the shock-layer region from 4 to 4+ 1 - For k = 1 to AT 

L Form the initial condition in newly refined regions. 

2. Use linear interpolation to compute the ghost values of Wq 

3. March Wo to based on scheme (3.7), 

IV. Project W 0 onto W 0 . 

V. Correct values Wq L and Wq r based on the shock-layer fluxes. 

Algorithm 1 Numerical Method . 


side of the refined region. There are a few points that need to be clarified in this algorithm. 
The interpolation to obtain ghost values (i.e. W^) is bi-linear interpolation based on 1 , 

Wql-i and Wq The initial condition for this problem is derived by imposing mass 

conservation; thus, the fine-grid values are all initialized to the value of the solution at the 
cell center. Improvements in the initialization procedure is a subject of further research. 

The correction of the coarse-grid values in Step VI is to use the same discretization that 
was used when the values were originally computed, but to modify the fluxes at the boundary 
of the domain to reflect what happened on the refined region. That is, to update Wq we 
would use scheme (3.5) with (3.6) for but we would compute with the formula 


F k ~ x , 

1 L+1/2 


1 


K - 1 




H- 


One may verify that this results in a globally conservative method. Also, this treatment 
of the boundary is consistent with the boundary conditions imposed in §2.1. Namely, this 
treatment of the interface is consistent with the viscous equation from problem (2.1) being 
the model at the interface between the subdomains. 


4.2 Dirichlet Treatment of Interface 

As a comparison to the flux boundary condition, we also implemented the heterogeneous 
domain decomposition method with dirichlet boundary conditions at the interface. This 
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is an interesting comparison, since there was little difference in the results when the two 
different treatments of the interface were used (this is discussed in §6). 


5 Detection of Interface 


We present the detection of the interface for the sake of completeness. Detection of the 
interface based on computational data results in a method that can have a different location 
of the internal-layer subdomain for each time step. The detection for the numerical method 
is based on obtaining an approximation to 


dWo dF(W 0 ) d_ 

dt ^ dx dx 



This term is the residual from using Wo as an approximation to the solution of (2.1). The 
residual is of magnitude (^(Ax" 1 ) in either a shock layer or in a zone where a shock interacts 
with other singularities. 

It is also possible to use an approximation of the viscous term j^Wq^^k) to localize 
some of the singularities. For example, this viscous term will be of order (^(Axr 1 ) in a 
shock layer or in a zone of interaction. This method is not as reliable as using the residual, 
however. Other types of behavior can be located and identified using these techniques [7]. 


0 Application to the Isentropic Gasdynamic Equa- 
tions 

In this section we examine the interface treatments on the viscous isentropic gasdynamic 
equations 

dudv 
dt dx 

dv d / 1 \ _ d ( 
dt dxxu^J € dx \dx ) ' 

Here u is the inverse of the density and v is the velocity. These equations are obtained from 
the conservation of mass and momentum in Lagrangian coordinates assuming that u is equal 
to the pressure raised to the power (the perfect gas law). The experiments were run 

with 7 = 2.2. 

The problem is a right-traveling shock interacting with and a left-traveling rarefaction, 
both of which eminate from the origin. An analytic self-similar solution (a rarefaction emi- 
nating from the origin) to the inviscid isentropic gasdynamic equations is given by 

= m 
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+ const . 


_ 2 7 1 /( 7 + 1 ) 

*■(*'') = -7TT-(j) 


(6.9) 


An initial condition with a shock and rarefaction eminating from the origin is constructed by 
connecting left values to middle values with a rarefaction. The middle values are connected 
to the right values with a shock. Thus, the initial condition is given by 

Ul> for x < 0 
Ur, for x > 0 


u(x,0) = | 
v(x, 0) = 


Vl, for x < 0 
Vr, for x > 0 


( 6 . 10 ) 

( 6 . 11 ) 


where 


U L = 1.4709, U R = 2.5000, V L = 1.0388, V R = 0.8050. 


The middle value of the solution between the shock and rarefaction is (17 m, Vm) = (1.973, 1.356) 
We remark that the middle values were was chosen using the Rankine-Hugoniot condition 

Vm - Vn l /Ur - 1 JIP M 

Ur-Um V r -V m • 

We expect the the viscous perturbation to have little or no effect on the speed at which 
shocks and rarefactions travel; thus, we will compare the viscous solutions to the solutions 
given above. 

The method was run with t — .01. The discretization parameters for numerical solution 
in the outer region have CFL number At/ Ax = .1, and Ax = .02. The discretization on 
the scaled coordinates inside the shock-layer is based on A£ = .1, with the CFL condition 
Ar/A£ < .025 and the stability condition Ar/A£ 2 < .1. These values are well within the 
limits imposed for the stability of the finite difference methods. 

Figure 6.2 depicts the evolution of the internal-layer subdomain when the two differ- 
ent boundary conditions are used. The errors generated by using the dirichlet boundary 
condition when the rarefaction is trying to exit the internal-layer subdomain result in a 
larger computed second derivative, and the detection scheme kept the rarefaction inside the 
internal-layer much longer. The solution projected onto the coarse grid at the end of the 
computations showed little difference between the two methods (Fig. 6.3). The primary dif- 
ference is the visual artificats at the boundary of the internal-layer subdomain at the point 
when the rarefaction is exiting the subdomain (Fig. 6.4). 


7 Conclusion 

Clearly the best interface condition is the flux-based treatment; however, the dirichlet bound- 
ary conditions did not induce as many errors as expected. One explaination of the lack of 
errors may be that the internal-layer subdomain boundary moves fast enough that waves 
propagating out of the internal-layer subdomain are allowed to pass across the bbundary 
by the oscillations in the boundary. More studies are planned with the goal to identify the 
precise nature of the errors associated with the interface treatments. 
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Time Step 


Evolution of IL (Dirichlet) Evolution of IL (Flux) 



- ! .0 -0,5 0.0 Oft 10 -l.n -0.5 0.0 05 10 

Spatial Location Spatial Location 


Figure 6.2: Evolution of the Internal-layer Subdomain. 


Solutions at t = .24 (Flux) Solutions at t = .24 (Dirichlet) 



Figure 6.3: Solution on Coarse Mesh at t — .24. 
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Solutions at t = .16 (Flux) 



Solutions at t = .16 (Dirichlet) 



Figure 6.4: Solution on Fine Mesh at t — .16. 
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